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Abstract. We study the inverse problem of determining the residual stress in Man's 
model using tomographic data. Theoretically, the tomographic data are obtained at the 
zero approximation of geometrical optics for Man's residual stress model. For compres- 
sional waves, the inverse problem is equivalent to the problem of inverting the longitu- 
dinal ray transform of a symmetric tensor field. For shear waves, the inverse problem, 
after the linearization, leads to another integral geometry operator which is called the 
>s mixed ray transform. Under some restrictions on coefficients, we are able to prove the 

uniqueness results in these two cases. 



1. Introduction 

Residual stresses are the stresses existing within an elastic body in the absence of 
external loads. They can be caused, for example, by heat treatment. The (Cauchy) 
residual stress tensor R is symmetric, Rij = Rji, satisfies the equilibrium equation 

div^ = in M (1.1) 

and the traction-free boundary condition 

RN = on dM, (1.2) 

where M C K 3 is the domain occupied by an elastic body and N is the normal to 
the boundary dM. Uncontrolled residual stresses are often detrimental to the life of 
a structural component. Nevertheless, in some cases, residual stresses are introduced on 
purpose to make materials or structures more resistant to damage, for example, toughened 
glass. It is therefore desirable to design reliable methods to determine the residual stress 
in a body. Besides its practical value, the determination of residual stresses gives rise 
to many challenging mathematical questions as well. Based on different model equations 
for the residual stress and using different types of measurements, several versions of the 
inverse problem of unique determination of the residual stress were studied in [HU] . [Ho] . 
[IWYlj . |IWY2j . |IMN] . [Ha] (see also [Rol], [Ro2] for related results). 

In this paper, we study the problem of recovering the residual stress tensor from the 
viewpoint of tensor tomography. To this end, we consider the time-stationary elastic 
system with residual stress that is assumed to be comparable with the inverse of the 
frequency. To be more precise, we set R = R/co, where welis the (angular) frequency 
of an elastic wave. In what follows, we refer to R as the residual stress tensor again 
although its physical dimension is time-stress = time-force/area. The time-stationary 
elastic wave equation is 

diva + co 2 pu = 0, (1.3) 
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where u is the displacement vector, p > is the density and a is the first Piola-Kirchhoff 
stress tensor. In view of Man's model |Maj . the first Piola-Kirchhoff stress tensor is 
expressed through the displacement vector by the equation 



a = A(tr£)J+ 2fie + - 



R + Vu-R 

(1.4) 

+ ^(tr e)(tr R) I + v 2 (tr R)e + v 3 ((ti e)R + (treR)I) +u 4 (eR + Re) , 

where e is the infinitesimal strain tensor, A and \x are Lame's parameters, and fi, . . . , V4 
are Man's parameters. Ignoring the term with 1/oj in (11.41) . we obtain the standard 
Hooke's law 

a = A(tre)/ + 2/x£. (1.5) 

We refer to (11.51) as the background isotropic medium for (ll.4p . while fll.4jl is called the 
quasi-isotropic perturbation of the medium (11.51) caused by the residual stress tensor R. 
Here, the term "quasi-isotropic" is used to emphasize that the perturbation smallness 
parameter 1/oj coincides with the reciprocal of the frequency. The concept of the quasi- 
isotropic approximation or quasi-isotropic media was first introduced by Kravtsov |Krj 
for the Maxwell equations (see also |KO] ). This approximation is based on the method 
of geometrical optics with rays generated by the background isotropic medium. Its gen- 
eralization to elastic waves is presented in \Sh\ Ch. 7]. 

Similar to the method used in [Sh] for a general anisotropic perturbation of an isotropic 
elastic medium, we start with applying the classical ray method to (11.31) with the stress 
tensor given by (11.41) . We restrict ourselves to considering the zero approximation of 
geometrical optics. As compared with the classical case of isotropic media, our formulas 
for the zero approximation have two distinct features. First, an additional factor appears 
in the formula for the amplitude of a compressional wave to describe the accumulation of 
the wave phase along a ray due to the residual stress. Second, the Rytov law for shear 
waves has got an additional term depending linearly on the residual stress. The inverse 
problem we study in this work is to determine the residual stress from the results of 
registration of the compressional or shear waves on the boundary of the medium under 
investigation. For compressional waves, the inverse problem is equivalent to the problem 
of inverting the longitudinal ray transform of a symmetric tensor field. For shear waves, 
the inverse problem, after the linearization, leads to another integral geometry operator 
that is called the mixed ray transform. 

The paper is organized as follows. In Section 2, we discuss the quasi-isotropic approxi- 
mation for Man's model of residual stress. For compressional waves, we derive the formula 
for the amplitude; while for shear waves, we obtain the Rytov law. All the content of 
Section 2 is actually some modification of the corresponding arguments from \Sh\ Section 
7.1]. We have chosen the following compromise presentation style in Section 2: all details 
are presented for compressional waves while our arguments for shear waves are condensed 
as far as possible. In Section 3, we investigate the inverse problem of determining the 
residual stress from measurements of compressional waves on the boundary. The inverse 
problem of determining the residual stress using shear waves is discussed in Section 4. 

2. QUASI-ISOTROPIC APPROXIMATION FOR RESIDUAL STRESSES 

In this section, we derive the quasi-isotropic approximation for the residual stress model 
1.41) in detail. We first write (11.41) in the form 



a = A(tr e)I + 2fie + ~(R + Vu ■ R + ce), (2.1) 

cj 
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where c = c(R) is a fourth rank tensor depending linearly on R. To write down the 
dependence explicitly, we reproduce (jl.4p and (12. 1 p in Cartesian coordinates 



Ojk — ^ £ ppfijk + 2/i£jfc H f-Rjfc + VjWp • R p k 

-\- V\Rpp£ qq& jk VlRpp^- jk VsiEppRjk ~\- £pqRqp$jk) ~\~ f£\E jpRpk Rjp^-pk) 



(2.2) 



and 

Cjk = teppSjk + 2fJ>£jk H — (-Rjfc + Vjtip • -Rpfc + Cjki m £ im)- (2.3) 
Comparing these two formulas, we deduce 



Cjkim = v\Rppbjkbi m + —Rpp{bji&km + 8j m 8ki) 

+l' 3 (Rj k 5i m + Rlmfijk) + — (Rjlfikm + Rjmfikl + Rkl^jm + Rkm^jl)- 



(2.4) 



Since we are going to use curvilinear coordinates, we rewrite (12. 4p and (I2.3P in the 
covariant form 

Cjklm ^1 

(tr E R)gjkgim + —(tr E R)(gjigkm + gj m gki) 

v, (2-5) 

+ y^{Rjkgim + Rlm9jk) + -^(RjWkm + Rjmgkl + RkWjm + Rkm9jl), 

a jk = \e p g jk + 2fie jk + -{R jk + u p .jR pk + c jk i m e lm ). (2.6) 

Here (gj k ) is the Euclidean metric tensor such that \dx\ 2 E = gj k dxidx k , (tr E R) = gi k Rj k) 
£ im _ g ji g km £ , k ^ and £ m _ g im £im _ Recall that e im = (u l;m + u m .j)/2, u x = gi k u h , 

dui k ; du l , k 

u h m = ^— - - l-im u k and u m = — — + r fcm tt , 



n v/-\ 1- r\ I hrinr a 



where T l km are the Christoffel symbol given by 

km 2 g y dxm t- dxk dxj )■ 

In curvilinear coordinates, ( II. 3p is expressed as 

a jk . k + uj 2 pUj = 0, (2.7) 

where a jk . = g lm a jk . m and 

jk; m ~ Ihf" ~ jm ~ km j ' 
It is clear that the tensor c possesses the symmetries 

Cjklm C-kjlm C-jkml Qrnj'fc (2-8) 

as follows from ( 12. 5p . The equilibrium equation ( II. ip is now written as 

R jk - k = 0. (2.9) 

We consider propagation of small elastic waves in a medium described by ( 12. 7p with 
the constitutive law (12. 6p . We reproduce some arguments of book |Sh] . The method of 
geometrical optics consists of representing a solution to the system by the asymptotic 
series 

oo m oo m oo m 

iu)T \^ M J _ iwr \^ £ jk __ iuir \^ ^jfc 



£ jk~e 2^ (4.„\m> a Jk~e ^ 



U ' = 6 

m=0 v ' m=— 1 v 
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where r = t(x) is a real function (eikonal). We insert the series into the equations under 
consideration, implement differentiations and equate the coefficients at the same powers 
of the frequency u on the left- and right-hand sides of the so-obtained equalities. In such 
a way we arrive at the infinite system of equations 

TYL L f 7TI 77% 7Tl~\~ 1 TTh I 1 \ 

£ jk = ^[ U j;k + U k . j + U j T- k + U k T.jj (to = -1, 0, • • • ), (2.10) 

m h m+1 k m+2 . , . . 

cr jk ; + VjkT- k -p Uj = (m= -2,-1,---), (2.11) 

a ife = A e^-fc + 2fie jk , (2.12) 

a jk = \e p p g jk + 2pe jk + i(T ;j R kp u p + R kp u ^ + c jkpq e pq ) (m = 0, 2, 3, ...), (2.13) 

a jk = Xe p g jk + 2pe jk + i (R jk + r. jR kp u p + Rkjfcj + c jkpq e pq ), (2. 14) 

— 1 —2 —2 —1 

where it is assumed that u = e = a = 0. Observe that a is a symmetric tensor while 
er (m > 0) is not symmetric. Putting m = — 1 in (I2.10p and m = —2 in (12.111) . we have 

-i 1 /Q o \ . , 

5 ifc= 2 V jT; fc + MfcT; V ( ^ 

and 

4^-/3^=0. (2.16) 

Observe that the residual stress does not participate in (I2.12p . (12.151) . and (12.161) . This 
means that geometry of rays is the same for the background isotropic medium (ll.5p and 
quasi-isotropic medium (ll.4p . As well known (see, for example, |Sh| p. 702]), (12.121) and 
(I2.15p - (l2.16p imply that r satisfies the eikonal equation 

|2 2 



where 



Vr 



2 2 P 2 2 P 

n = n n = or n = n= —. 

p A + 2/i jj, 



Next, we calculate the (complex) amplitude A p of a compressional wave by repeating 
arguments of [Shi Section 7.1.3]. Recall that the amplitude A p is defined by 

o A + 2/i o _ . _ . Vr 
u = — (u, Vr)Vr = A p - 



P v ' ' HVr| 

We fix a solution r to the eikonal equation with n = n p and introduce ray coordinates 
in a neighborhood of a ray, i.e., curvilinear coordinates 3 such that x 3 = t and 

the coordinates surfaces orthogonal to the coordinates lines X — Xqj X — 

that are geodesies of the metric ds 2 = n 2 gj k dx j dx k . In such coordinates 

a a p = \ng a pA p , a a3 = 0, a 33 = (A + 2p)n~ 1 A p (2.17) 

(see jSh, formulas (7.1.23)]). Here Greek indices assume the values 1,2. Likewise, we can 
get that 

o 1 /o o \ o o 1 , nin , 

£ a p= ^ y u or,p + U P;a) > ^33=%; 3 + %, (2. 18) 

a 33 = A^V^ + (A + 2pi)e 33 + i(ni2 33 A p + c 33p(? l: lw ) (2.19) 
*s* ; " + n 2 a 33 -p l 3 = 0. (2.20) 
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Formulas (7. 1.27)- (7. 1.29) of |Sh] remain unchanged, i.e. 



-i 



1 dg a p o ^dAp 
rn—^Ap, u 3 . 3 = n 



2 dr 

dg a fs 
dr 



fin-^Ap, 



^ 3;3 = (A + 2 / iK 1 ^4 



dr 



Ot 



i d(X + 2a) , 9 <9n 



(2.21) 



.4 



0"3fc: 



ax 



(A + 2jj)n-Q— + [ jing 



a adg a p , d{\ + 2p) 



dr 



+ n 



<9r 



On \ 

+ (X + 2fi)—jA pj 



o 1 dg a p o _!<9v4 P , i 



Using fLT23|) in fl2TT9|) implies 

^33 = (A + 2\i)n 



dr 



(2.22) 
(2.23) 



i 



' -i dA P + 1 XnT^^A, + (A + 2/i)w 3 + i(«i?Bs^ + c 33p(? e 



<9r 



Or 



Inserting this value for cr 33 and value (12. 22ft for a 3 j. . k into (I2.20p . we arrive at the relation 



Ft A 

2(A + 2/i)n^ + 

OT 



-(A + 2[i)ng 



dr 



+ n 



Or 



+ (A + 2^' 



A r 



+ [(A + 2p)n 2 - p}u 3 + i(n 3 R 33 A p + n 2 c 33pq e pq ) = 0. 



In view of n 2 = p/(A + 2/i), the coefficient at u 3 in this formula is equal to zero. Thus, 
inserting into the last formula the expressions 



-J-a/3 



L a3 



0. 



n 3 A r 



that follows from (7.1.11), (5.1.20) and (7.1.22) of |Shj . we arrive at the equation for the 
amplitude A p 



+ 



dA p ( g aP dg a p 1 <9(A+2/i) 1 dn nr 



+ 



+ 



+ 



#33 + 



in 



dr ' V 4 dr 2(A + 2/i) dr 2n dr X + 2p 2(X + 2p) 

Using notation (5.1.22) and formula (5.1.23) of [Shj . we write this in the form 

d 



C3333 ) A p = 0. 



rr 



R33 ~ i 



n 



2(A + 2p) 



c 3333) 



which, together with n 2 = n 2 , gives 

d_ 



In [A p ^ Jpvp 



1 1 

#33 — 7. r C 3333, 



2pv% 



where J is the geometrical spreading (see (Shi (5.1.22)]) and v p = l/n p is the velocity of 
compressional waves. This implies the following analogous of [Sh| formula (7.1.32)]: 



.4 



C 



\/ j P v p 



exp 



1 

PVp 



1 



2pv% 



y^ k 4^r dr 



gCjklmT'y 7 7 



(2.24) 



where 7 J = c?7 J /dr and C is a constant for a given ray 7. 

Next, we express the second integrand through R. Since 7(7-) is a geodesic of the 
Riemannian metric \dx\ 2 = v~ 2 \dx\ 2 E) it satisfies 



9jki 3 i k = |7ll = v. 



p I'ylp 
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under the assumption that 7 is parameterized by the r-length, i.e., \j\ 2 = 1. Taking this 
into account, we obtain from (I2.5P 

c jk imi j i k i l i m = v 2 p (2(1/3 + ^)R jk i j i k + vfa + ^)tr B r) . (2.25) 

Here trg R is the trace of R with respect to the Euclidean metric (gjk), i-e., tr# R = gi k Rj k . 
Since 7 is a geodesic of the Riemannian metric hj k = v~ 2 gj k , it is more natural to use the 
Riemannian trace tr R = h? Rj k = fptrgi?. Thus, (I2.25p takes the form 

CjMm'f^T = v 2 p (2(1/3 + ^)Rjki j i k + {v x + U 2 )tT i?) . (2.26) 

Substituting this value into (I2.24p . we obtain the final formula for the amplitude of the 
compressional wave 



C_ 

\/Jpv P 



A p = — - = exp 



-ij ^4 ((1 + ^3 + ^)R 3 ki j i k + \{vx + u 2 )tiRj dr 



(2.27) 



The physical meaning of the formula is as follows: the residual stress R distorts the phase 
of the compressional wave. 

Finally, we consider propagation of shear waves by repeating arguments of Sections 
7.1.4 and 7.1.5 of |Sh]. It can be easily checked that the residual stress R does not change 
equations (7.1.43) of |Sh] . Therefore the same formula A s = CjsJ Jpv s is valid for the 
amplitude of the shear wave as for the background isotropic medium where v s = l/n s 
is the velocity of shear waves. Moreover, Rytov's law in |Sh| formula (7.1.51)] does not 
change either. We reproduce the formula here 



dr J . "4pv 



-/-^ - i 3 i q )c qklm i k i m r]\ (2.28) 







where 77 is the polarization vector defined by Uj = A s n~ 1 r]j. Here D/dr = 7 J 'Vj with 
the covariant derivative taken with respect to the Riemannian metric dr 2 = v~ 2 \dx\ 2 E = 
hjkdx J dx k and indices are raised with the help of the same metric, i.e., 7 fc = h^ kA fj and 
f] 1 = y l 7]j. 

3. The inverse problem for compressional waves 

First of all we emphasize the following feature of our approach. While considering the 
forward problem, we impose no boundary condition on the displacement vector u at the 
boundary of the domain under consideration. Thus, we treat the problems as if the waves 
propagate in an unbounded medium, and use the boundary only as a surface at which the 
sources and detectors of oscillations are disposed. In fact, due to the reflection effects on 
the boundary, the possibility of registration of information that is used below as the data 
for inverse problems seems to be rather problematic. Here we will not settle this question 
but only attract reader's attention to the fact of its existence. 

In contrast to the content of the last paragraph, the following remark gives the possibil- 
ity of measuring the data for compressional waves regardless to any boundary condition. 
There exists an alternative version of the geometrical optics method which is based on 
the analysis of propagation of the wave front of a non-stationary elastic wave, see |CJj . 
In this version, the integral participating on (I2.27P appears as the first order perturbation 
for the travel time of a compressional wave. More precisely, for a fixed geodesic 7 of the 
metric \dx\ 2 = v~ 2 \dx\ 2 E between two boundary points, let Ti/ W (7) be the propagation 
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time of the compressional wave along 7 in the quasi-isotropic medium (11 .4p and To (7) be 
the corresponding travel time in the background isotropic medium (jl.5p . Then 

Ti/M-T ( 7 ) = ~D(i) + o(^), 

where 

D (l) = [ A ((! + ^ + ^) R jki j l k + 75(^1 + nOtrfl) dr (3.1) 

7 

is just the integral participating in (I2.27p . The corresponding result is obtained in |CJj 
in the case of a perturbation of the form 

0~jk — ^£pp3jk + 2/i£jfc H Cjklm^lm- 

UJ 

By repeating arguments of [CJJ, one easily sees that the result is true if the last formula 
is replaced with (12. 3p . Thus, data (13. ip can be obtained by measuring travel times for 
compressional waves. Let us consider an elastic wave initiated by a 5-kind source at the 
initial point of 7 which starts at the time t — 0. The wave will be a mixture of different 
body and surface waves including secondary waves caused by reflections at the boundary. 
Nevertheless, Ti/ W (7) is the first arrival time to the final point of 7 since compressional 
waves are the fastest elastic waves. In the simplest case of constant parameters A, p, and 
p, To (7) is equal, up to a constant factor, to the Euclidean distance between the endpoints 
of the straightline segment 7. So, the only problem is the sufficiently precise measurement 
of the travel time ^^(7). 

Studying the inverse problem, we will first consider the case of constant coefficients 
A, p,p, Vi, ... ,1/4 since this case is much easier than the general one and, most probably, 
this case is of the most importance for applications. 

For the inverse problems studied here, we assume the material parameters X,p,p, 
ui, ... ,1/4 to be given a priori and only the residual stress R to be unknown. In prac- 
tice, some of these parameters are also unknowns to be determined. Therefore the inverse 
problem of recovering residual stresses, when some of material parameters A, pi, p, V\, . . . , V4 
are also unknowns to be recovered, is also worth of investigation. But this new inverse 
problem is much harder because it is a nonlinear problem. 

3.1. The case of constant coefficients. Let the medium under consideration be con- 
tained in a bounded convex domain McK 3 with smooth boundary dM and let each of 
the material parameters A, p, p, V\, . . . , z/ 4 be constant. In this case, the Riemannian met- 
ric h = v~ 2 \dx\ 2 E is a constant multiple of the Euclidean metric g = \dx\ E and geodesies 
are intersections of straight lines with M. The equilibrium equation (12 .9p means that R 
is a solenoidal tensor field. We extend R to the whole of M 3 by zero outside M. Then 
the extended tensor field R is solenoidal on the whole of M. 3 in virtue of the boundary 
condition (II. 2p . 

We study the inverse problem of recovering the residual stress tensor field. To this 
end, assume that we can dispose a source of compressional waves at every point of the 
boundary dM and measure the phase of a compressional wave on the same surface dM. In 
virtue of (I2.27p . our data are integrals (13 .ip that are known for every line 7 of M 3 . Initially 
in (13. ip . 7 is parameterized by the arc length in the metric h, i.e., \ j\E = v p = const and 
the trace is understood with respect to the metric h, i.e., tr R = v 2 tr E R. After a simple 
rescaling, we obtain the same formula (13. ip . where now I'jIe = 1 and tri? is replaced by 
the Euclidean trace tr eR- From now on in this subsection, we use the Euclidean metric 
only and write tr R instead of trg-R. 
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To get a well defined inverse problem, we have to impose some restrictions on the 
material parameters. Indeed, if for example V\ + v<i = —2(1 + z/3 + ^4)/3, then the 
integrand on (13. ip is identically equal to zero for R = g. Therefore, in this section, we 
assume that 

3(i/i + u 2 ) + 2(1 + u 3 + 1/4) + 0. (3.2) 

The first term of the integrand on (13.11) is considered as the leading term. Therefore, in 
this section, we also assume that 

1 + 1/3 + 1/4 ^ 0. (3.3) 

Introducing the notations 

1 + 1/3 + 1/4 z/i + i/ 2 . . 

t = a R, a = — ; r 3.4 

J pv* ' 2(l + z/ 3 + z/ 4 ) 1 ; 

(/ is zero outside M), we write data ( 13.1 ft as 

(/(/ + a(tr/) ))( 7 ) = / (f lk + a(tif)g lk )^fdr, (3.5) 



where I is the (longitudinal) ray transform on IR 3 which is defined in Section 2.1 of |Shj . 

Let us remind the theorem on decomposition of a tensor field into solenoidal and po- 
tential parts (Theorem 2.6.3 of |Sh] ) : every symmetric tensor field u = (ujk) € L 2 on R 3 
can be uniquely represented as 



u = u + dv, 5u = 



where v is continuous outside M and satisfies v(x) — > as \x\ — > 00. Here the operators d 
(inner derivative) and 5 (divergence) are defined in local coordinates by formulas (dv)jk = 
(VjVk + VkVj)/2 and (5u) j = g kl VkUji respectively, V being the covariant derivative with 
respect to the Euclidean metric g. The summands u and dv of the decomposition are 
called the solenoidal and potential parts of the tensor field u respectively. 

We now investigate the question of uniqueness of a solution to the inverse problem. Let 
a solenoidal field / satisfy /(/ + a(tr f)g) = 0, a = const. By Theorem 2.15.1 of [Sh], this 
means that the solenoidal part of / + a(tr f)g is equal to zero. Since the solenoidal part 
of / coincides with /, we obtain the equation 

f + aS((trf)g) = 0, (3.6) 

where S(u) denotes the solenoidal part of a tensor field u. Recall that the Fourier trans- 
form interweaves the operators S and T, where Tu stands for the tangential part of u, 
see Section 2.6 of |Sh] for details. Applying the Fourier transform to the last equation, 
we obtain 

/ + a(tr/> = 0, (3.7) 

where the tensor field e — T(g) is expressed in Cartesian coordinates by £jk(y) = 8jk — 
VjUk/lyl 2 , V being the variable in the Fourier space. 

Applying the operator tr to equation (13.71) . we obtain (1 + 2a)tr / = 0. Thus, the 
inequality 1 + 2a 7^ is the necessary and sufficient condition for the uniqueness of a 
solution to the inverse problem. Recalling (13. 4p . we write the condition as 

V\ + Z/2 + ^3 + ^4 7^ "I- (3.8) 

If (13. 8p holds, equation (1 + 2a)tr/ = gives tr / = 0. Then (13. 7p implies / = and 
therefore / = 0. 
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Under hypothesis (|3.8|) . an explicit inversion formula for recovering a solenoidal tensor 
field / from the data 7(/ + a(tr f)g) can be easily derived from the corresponding inversion 
formula for /, see Theorem 2.12.2 of |Shj . The corresponding stability estimate can be 
also obtained on the base of the Plancherel formula for the ray transform, see Section 
2.15 of |Shj . Moreover, to recover a solenoidal /, we do not need to measure ray integrals 
(I(f + a(tr f)g)) (7) for all lines 7 of M 3 . Repeating arguments of |Sh2] . we see that three 
families of lines are sufficient for an effective reconstruction algorithm, each family consists 
of all lines parallel to a coordinate plane. 

If ( 13. 8 p does not hold, i.e., if v\ + v?, + + v± = — 1, then the space of solenoidal tensor 
fields / satisfying !(/ + a(tr f)g) = can be explicitly described. Indeed, in this case tr / 
can be an arbitrary function and equation ( 13. 6 j) gives / = S(ag) with an arbitrary scalar 
function a. 

3.2. The case of variable coefficients. Let again M C R 3 be a closed bounded domain 
with smooth boundary dM. Now the material parameters A, p, p, U\, . . . , 1/4 are assumed 
to be known smooth functions of a point x G M. Let v p = a/(A + 2/x)/p be the velocity 
of compressional waves. By g we denote the Euclidean metric and by h = v~ 2 g, the 
Riemannian metric corresponding to compressional waves. Assume (M, h) to be a convex 
non-trapping manifold (CNTM) in the sense of definition given in |Shlj . The same defi- 
nition is presented in Section 4.1 of [Sh] but the term "compact dissipative Riemannian 
manifold" is used instead of CNTM. 

We consider the inverse problem of recovering the residual stress tensor field. Our data 
are integrals (13.11) that are known for every geodesic 7 of the metric h with endpoints in 
dM. The geodesic is parameterized by the arc length in the metric h, i.e., I7I = 1. By the 
same arguments as above, we assume inequalities (13.21) and (I3.3P to be valid everywhere 
in M. On using the same notations (13.41) . we write data (13.11) as 

/(/ + a(tr f)h) = J (f jk + a(tr f)h jk )f^ k dr, (3.9) 

7 

where I is the (longitudinal) ray transform on the CNTM (M, h) which is defined in 
Section 4.2 of [SE]. 

For a compact Riemannian manifold (M, h), the theorem on decomposition of a tensor 
field into solenoidal and potential parts (Theorem 3.3.2 of |Shj ) is valid in the following 
form: every symmetric tensor field u = (ujk) can be uniquely represented as 

u = u + dv, 5u = 0, v\qm — 

where the operators d (inner derivative) and 5 (divergence) are defined in local coordinates 
by the same formulas {dv)jk = (V/Wfc+Vfc'Uj)/2 and (Su)j = h kl VkUji respectively, V being 
the covariant derivative with respect to the metric h. 

According to Theorem 4.3.3 of [Sh] . the solenoidal part of the tensor field / + a(tr f)h 
can be uniquely recovered from data (I3.9P under the assumption 

k + (M,h) < 1/3, (3.10) 

where k + (M,h) is some curvature characteristic of the CNTM (M,h). For the metric 
h = v~ 2 g, condition (I3.10P holds if the function v p is sufficiently C 2 -close to a constant, 
the degree of the closeness depends on the size of the domain M. For such a manifold, 
the null-space of / consists of potential fields that can be represented in the form dv with 
a covector field v satisfying the boundary condition v\sm = 0. 
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For variable coefficients, the main difficulty relates to the following circumstance: the 
equilibrium condition (I2.9P does not mean anymore that / is a solenoidal tensor field with 
respect to the metric h. Indeed, (I2.9P can be rewritten in terms of / as 

4 

6 E (bf) = with b= ^ , (3.11) 

^ ' 1 + ^3 + ^4 

where 5e is the divergence with respect to the Euclidean metric g. We are going to prove, 
at least under some restrictions on the material parameters A, /i, p, v%, . . . , 1/4, that a tensor 
field / is uniquely determined by data (13.91) if it satisfies (13. lip . 

First of all we will rewrite equation (13.1 ip in terms of the divergence Sf with respect to 
the metric h. Denote c = v%, then g = ch. On using standard formulas of tensor analysis, 
one easily calculates 

= bc-\5f) j + c- 1 (f jk V k b- ^c- 1 ((tr/)V J c-/ jfc V fe c)), 

where the covariant derivative and trace are understood with respect to the metric h. 
Therefore (13. lip is equivalent to the equation 

(5/), + 6- x /ifcV fc 6 + ^trVifcVc - I c -i(tr f)V jC = 0. 

Denoting 

a = b~ l Vb + ^c _1 Vc, p = --c^Vc, (3.12) 

we write the equation in the coordinate free form 

6f + fa + (trf)f3 = 0. (3.13) 

We now investigate the uniqueness question to the inverse problem. Let a tensor field 
/ be such that I(f + a(tr f)h) = 0. We assume that the Riemannian manifold (M, h) is 
a CNTM and satisfies the curvature condition 1 13. 10p . Then, by Theorem 4.3.3 of [Shj . 
/ + a(tr f)h must be a potential field, i.e., there exists a covector field v on M satisfying 
the boundary condition v\qm = such that 

/ + a(tr f)h = dv. 

Taking the trace of both parts, we obtain 

(l + 3a)(tr/) = 5v. 

By (13. 2p . 1 + 3a does not vanish on M and we can write 

tr / = (1 + 3a)' 1 5v, f = dv - a(l + 3a)-\5v)h. 

Substituting these values into ( 13.131) . we arrive at the boundary value problem on the 
covector field v 

[ <TT3^ " M)V - + TTTa {S < aa ~ ' + (1 + ^TTTa') = °' (3.14) 

[ v\dM = 0. 

We have thus proved 

Proposition 3.1. Let a three-dimensional CNTM (M, h) satisfy the curvature condition 
A3.10\) and let a, a, and (3 be defined by l[3.4\ ) and A3.12\) . Every symmetric tensor field 
f satisfying A3.13\) can be uniquely recovered from data Ii3.9\) if and only if the boundary 
value problem flff. 1$ has no nontrivial solution. 
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Equation f l3.14|) is rather complicated in the case of general coefficients a, a, and 0. We 
are going to investigate the boundary value problem in the case when the background 
medium is sufficiently close to a homogeneous one, i.e., when a is close to a constant and 
a and /3 are small. Even in this case, we need to impose some restrictions on a. 

Let us first find a condition that guarantees ellipticity of the boundary value problem. 

Proposition 3.2. The operator jr^dd — 5d is elliptic if and only if the inequality 

|l + a|<|l + 3a| (3.15) 

holds on the whole of M. 

Proof. Recall that the principal symbols of the operators d and S are a/ — li% and \/ — lj% 
respectively, where i^ is the symmetric multiplication by the covector £ and j% is the 
contraction with £, see |Sht Section 3.3] for details. Therefore the principal symbol of our 
operator is 

By |Shl Lemma 3.3.3], j^i^ = |(|£| 2 i?+i^) on covectors, where E is the identity operator. 
Therefore 

°( T J ^ r d6-6d)=lm 2 E + ^) with k=±±±. (3.16) 
V 1 + 3a / 2 l + 3a 

On the other hand, for a covector v, 

m 2 E + K lac )v,v) = \^\ 2 \v\ 2 + K^,v) 2 . 

Hence, f)3. 15j) guarantees the positiveness of \£\ 2 E + Ki^j^ for all £ ^ 0. □ 

Under suitable assumptions on coefficients, the triviality of a solution to the boundary 
value problem (I3.14p is guaranteed by the following 

Theorem 3.3. Given a three-dimensional CNTM (M,h), let D be its diameter, i.e., the 
length of the longest geodesic. Assume the coefficients a, a, and (3 of equation (13. 14j) to 
satisfy 

3a + \a Q + ho + \{al + ffi)D 2 < 1, (3.17) 



2 2' 4 



where 

ao = sup 



ao = ^sup \a\J , (3q = ^sup 



aa — (5 



l + 3a 



1/2 



1/2 

r>n = I sun \n\ I 

l + 3a 

Then the boundary value problem (13 . 14j) has only trivial solution. 

The values of ao and /3q can be made arbitrary small by assuming the background 
medium to be sufficiently close to a homogeneous one. Therefore the main part of hy- 
pothesis ( I3.17|) is ao < 1/3 that is equivalent to 

a>-l/6 (3.18) 

It is interesting to compare (13.181) with the ellipticity condition (13. 15 j) . If— l/6<a<0, 
then the boundary value problem is not elliptic but has only trivial solution for sufficiently 
small a and /3. If both conditions ( 13. 15j) and (13.171) are satisfied, we can use the standard 
stability estimate for the elliptic boundary value problem with the trivial kernel to obtain 
some stability estimate in the inverse problem of recovering a tensor field / from data 
(13.91) . This stability estimate will be similar to that of |Sht Theorem 4.3.4]. We do not 
present it here. 

To prove Theorem 13.31 we need the following 
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Lemma 3.4. For a covector field v on a three-dimensional CNTM (M, h) satisfying 
v\dM = 0, the estimate 

\\v\\h < ^(2\\dv\\l 2 + \\5v\\l 2 ) 
holds where D is the diameter of (M, h) . 

Proof. Let QM = {(x, £) | x 6 M, £ e T X M } |£| 2 = = 1} be the unit sphere bundle. 

Introduce the L 2 -norm 

1 2 _ . / I../ M|2 



WAl^m) = J du x (£)dV(x), 

QM 

where du x is the volume form on the sphere Q x = QM C\T X M and dV is the Riemannian 
volume form on M. 

Given a covector field v on M, define two functions on QM 

<p(x, = Vi{x)C, ip{x, = Ulri.r)),^^. 

The functions are related by the equation Hip = ip where H is the differentiation with 
respect to the geodesic flow, see |Shl Section 4.4] for details. If v\qm = 0, then ip vanishes 
on the boundary of QM. Applying the Poincare inequality |Sh[ Lemma 4.5.1] with the 
weight A = 1, we obtain 

D 2 D 2 
y\\h(QM) < -yll^llx^nM) = -ylMli^nM)- (3.19) 

One can easily see that (compare with |Shl Lemma 7.4.2]) in the three-dimensional case 

II II 2 - II II 2 

W\\L 2 (nM) — g \\ V \\L 2 (M)i 

H\\h(nM) = ^WMhiM) + \\J dv \\h(M)) = Y^( 2 \\ dv \\h(M) + W^Wh(M))- 
Inserting these values into (13. 19ft . we get the statement of the lemma. □ 

Proof of Theorem 13.31 Let v be a solution to the boundary value problem (|3.14p . First 
of all we rewrite equation ( 13.141) in the form 

d( 5v) — Sdv — (dv)a + (6v)— — — = 0. 

Vl + 3a / K J K J l + 3a 

Take the L 2 -product of the equation with v, use the boundary condition v\qm = and 
the fact that d and —5 are dual operators (see \Sh\ Theorem 3.3.1]) 

IHI 2 = {j^-SvM + ((dv)a,v) - (^JLsv,v). (3.20) 
1 + 6a 1 + 6a 

We estimate each summand on the right-hand side of (I3.20p . Obviously 



a 



-5v, 8v) 



. ., <a \\5v\\ 2 . (3.21) 

1 + 6a 

The second summand on the right-hand side of (13.201) is estimated as follows: 

\((dv)a,v)\ < ( a ; /2 ||^||)(af || V ||) < ^a \\dv\\ 2 + a 3 \\v\\ 2 ). 
This implies with the help of Lemma 13.41 

\((dv)a,v)\ < ha °\f D2 \\dv\\ 2 + ^II^H 2 . (3.22) 
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The last term on the right-hand side of (13.201) is estimated in the same way: 
,aa — (3 

nr. i)r ) - 

10 " " 20 



-Sv, 5v) 



<^ W |. + jjft±£g| W |. (3.23) 



l + 3a 

With the help of (l3T2T]) - ([3723]) . f l3~20|) gives 

\\dvf < ^(5a + alD 2 + (3 3 D 2 )\\dv\\ 2 + (a + ^(10/3 + a*D 2 + $I> 2 )) HHP- (3.24) 

The inequality \\5v || 2 < 3||ck>|| 2 is obvious since 5v is the trace of the 3 x 3- matrix 
((dv)ij). With the help of the latter inequality (I3.24p gives 



1 - 3a - i« - ^A> - + Pl)D 2 



\dv\\ 2 < 0. 



The coefficient in the brackets is positive under hypothesis ( 13 . 1 T[) and the inequality 
implies dv = 0. With the help of the boundary condition v\qm = 0, this implies v = 0. □ 

4. The inverse problem for shear waves 



Let v s = \fjjjp be the velocity of shear waves. By g we denote the Euclidean metric 
and by h = v~ 2 g, the Riemannian metric corresponding to shear waves. Assume (M, h) to 
be a CNTM. We assume that, for every geodesic 7 : [0, 1] — > M, \j\ 2 = hj^j = 1 with 
endpoints in dM, we can activate the shear wave with an arbitrary initial polarization 
77(0) and can measure the final polarization r](l). Mathematically, this means that the 
fundamental matrix U (7) of system (12.281) is known for every geodesic 7 with endpoints 
in dM such that rj{l) = £7(7)77(0). We are going to recover the residual stress R from the 
data {7(7). 

Introduce the tensor field / by 

fjklm i ^pyQ \Cjlkm Cjmkl)- (4-1) 

It possesses the symmetries 

fjklm fkjlm fjkml flmjk 

as follows from (12. 81) . Express / through R by substituting ( 12. 51) into (14.11) 

fjklm = ~^ 4p t ,4 ( t/l (^ r R)(hjlhkm + ^jmhkl) + ~^ (t r R) ^^jkhlm + hjih km + hj m h kl ) 

+ l/ 3(Rjlhkm + Rjmhkl + Rklhjm) (4-2) 
+ -^-{2Rjkhl m + Rjlhkm + Rjmhkl + Rklhjm + Rkmhjl + 2Ri m hjk i 

Here tri? = W k Rj k is the Riemannian trace. 
Rytov's law (12.281) is written in terms of / as 

^)_ = (?-7if)W7V- (4-3) 

To write Rytov's law in a coordinate free form, we introduce the following notations. 
Define the linear operator fa by (far])j = fqrjk^Y^- F° r a nonzero vector £, let be 
the orthogonal projection onto ^ — {v \ (d 7 )) = = 0}- 111 coordinates (P^i])j = 
{$j ~ j^2^j^, p )Vp- Then the coordinate free form of ( 14.31) is 
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This equation can be also written in the form 



^ = (P.f^rj, 77(0) = ryo, (4.4) 



where, for a symmetric tensor u = (ujk), 

\w MP){sl ~w 

compare with |Sh| formula (5.2.1)]. is the orthogonal projection onto the subspace {u \ 
Ujk£, k = 0}. Let 7 _l (t) be the (complex) two-dimensional subspace of (the complexification 
of) the tangent space T 7(r )M consisting of vectors orthogonal to and let I® ,T : 

7 (0) — > 7 _l (t) be the parallel transport along 7 with respect to the metric h. The 
solution to the initial value problem (14. 4p can be written as 

r ] {l)=I^ l {U{ 1 )r ] o) for 77067^(0) (4.5) 

with some linear operator (the fundamental matrix) ^(7) : 7 _L (0) — > 7" L (0). Observe that 
[7(7) is a unitary operator since the matrix of system (I4.4j) is skew-Hermitian as is seen 
from (14. 2p . We consider the problem of recovering the tensor field R from the data U{^) 
known for all geodesies with endpoints in DM. 

The problem is strongly nonlinear since the data {7(7) depends on R in a nonlinear 
manner. Let us linearize the problem. To this end we represent the fundamental matrix 
by the Neumann series 

1 l T 

U( 1 ) = E + J F/{{P^){r))dT + j r/{{P^){r))dr j /f {{P^){t)) dt + . . . , 
00 
where E is the identity matrix, and delete the terms that are nonlinear in /. In other 
words, we use Born's approximation 

1 



U(i)-E= j ^°((P 7 / 7 )(r))rfr 





as the data for the linearized inverse problem. Since the integrand is a symmetric operator, 
the data are equivalent to the quadratic form 

1 

((17(7) - E) V (0), 77(0)) = j (((Pyfy(r)) V (r), V (r)) dr (4.6) 

'0 

on the two-dimensional vector space of vector fields r\ (r) that are orthogonal to 7 and 
parallel along 7 in the sense of the metric h, i.e., satisfy ^ = 0. We denote this space by 
7- 1 -. Since 

((Pjfa)v,V) = ( P i(fiV),v) = (fjV, p -rV) = (fjV,v) = fjkimi JA f k V l V m , 
(14.61) can be written as 

1 

((i7( 7 ) -E) V (0), 77(0)} = (Lf)( T ,77) := j /^ m (7(r))7 J (r)7 fc (r)r7 Z (r)r/ m (r) dr. (4.7) 



The operator L defined by this formula is called the mixed ray transform (compare with 
formula (7.1.55)]). 
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We express the integrand of (14. 7p through R by substituting value (14. 2[) for /. On using 
the relations |7| = 1 and (7,77) = 0, we obtain 

(Lf)(j,r)) = -1 J -^{v^rfrf + iy A R jk ^j k \r,\ 2 + z/ 2 (tr P)M 2 ) dr. (4.8) 



7 



The first term of the integrand on H4.8[) is considered as the leading term. Therefore we 
assume in this section that the function z/ 4 does not vanish in M. Note that we do not 
use the equilibrium condition (12 .9p in this section. Introducing the notations 

R, a = —, (4.9) 



4pv 4 s ' ~ z/ 4 
we write (14. 8p as 

(L/)(7,r/) = J F jk r 1 1 r 1 k dr+\r 1 \ 2 j F jk f j k dr + \rj\ 2 J a(txF)dr. (4.10) 

777 

The first term on the right-hand side of (14. lUp is the transverse ray transform ( JF){^f, 77), 
compare with |Sh| formula (5.1.72)]. The second term coincides, up to the factor |t7| 2 , with 
the longitudinal ray transform (IF)^), while the last term coincides with J(a(tr F)/i)(7). 
Therefore (I4.10p can be written as 

(JF)( 7 ,77) + M 2 (/F)( 7 ) + |r / | 2 /(a(trF)/ i )( 7 ) = (L/)( 7 ,ry) (77 G r 1 ). (4.11) 

We consider (14. lip as an equation in an unknown symmetric tensor field F = {Fj^.) while 
the right-hand side (Z//)(7,7?) is given. 

Equation (14. lip can be simplified. Indeed, let (771,772) be an orthonormal basis of 7- 1 -. 
Then 

(JF)(7,77i) + (JF)(7,r 72 )+2(/F)(7)+2J(a(trF)/ i )(7) = (L/)( 7 , r 7l ) + (L/)( 7 , 77 2 ). (4.12) 
Since 



(JF)( 7 ,r/i) + (JF) (7,773) = J (tr P^F)dr, 



(I4.12p can be written as 



I |2 f I |2 

' 7l htiP i F)dr+ |7?| 2 (LF)( 7 ) + \ri\ 2 1 (a(ti F) h) (7) 



(L/)(7,7 7l ) + (L/)( 7 ,7 ?2 ) . 



2 

Subtracting this equality from (14. lip , we obtain 

1 1 2 /> 1 1 2 

(JF)(7,7?) - ^- y (tr P^F)dr = (L/)( 7 ,t ? ) - PL ((L/)( T , 77a) + (L/)( 7l ^)). (4.13) 

7 

Let be the orthogonal projection of symmetric tensors onto the subspace {m = (ujk) \ 
Ujk£, k = 0, trw = 0} (compare with Section 6.2 of [Shj). The left-hand side of (I4.13P can 
be transformed as follows: 

I |2 r r I |2 

(JF)(7,t?)-^L J(trP 1 F)dT = J (((P^F)77, 77) - ^-(tr P^F)) dr 



7 



((Q^F)r),r)) dr=(KF)( 1 ,r ] ), 
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where K is the truncated transverse ray transform, see the definition in Section 6.2 of 
[Shj . Since (L/)(7, 77) is known for every 77 G 7 -1 , the right-hand side of (I4.13P is known 
too. Moreover, it is independent of the choice of an orthonormal basis (771,772) as follows 
from f)4.13p . Denoting the right-hand side of (14. 13j) by -D(7, 77), we arrive to the equation 

(KF) (7, 77) = £>( 7 , 77) (77 G r 1 ). (4.14) 

Let us distinguish the trace free part of the tensor F, i.e., represent it in the form 

F = F+ -(trF)h, where trF = 0. 
3 

By Theorem 6.6.2 of [Shj . the trace free tensor field F can be uniquely recovered from Kf 
if the CNTM (M, h) satisfies some curvature condition. Moreover, the stability estimate 

\\F\\ L 2<C\\KF\\ m 

holds with some constant C independent of F. 

On assuming F has been recovered, we can calculate IF and JF . Then 

(JF)( 7 ,77) + M 2 (/F)( 7 ) 

= ±(J((trF)h))( %V ) + M!(/((trF)/ i ))( 7 ) + (JF)( %V ) + |r/| 2 (/F)( 7 ) 

= ^(/((trF)/0)( 7 ) + iy^((tri^))(7) + (JF)(l,v) + \v\ 2 (IF)(l). 
Substituting this value into (14. lip , we obtain the equation 

(I((a + |)(trF))( 7 ) = j4(W)(7^) - (^)(7^)) - (/F)( 7 ). (4.15) 

The left-hand side of this equation is the ray transform of the scalar function (a + |)(tr F) 
while the right-hand side is known. By the way, the right-hand side must be independent 
of 1] G 7 as follows from the equation. 

By Mukhometov's theorem [Muj . the function (a + |)(trF) can be uniquely recovered 
from the ray transform J((a + |)(tr F)) if (M, h) is a simple manifold, see Section 1.1 of 
[Sh] for the definition of a simple manifold. Thus, the trace tr F can be recovered under 
the additional assumption that the function 3a + 2 does not vanish. Since the trace free 
part F has been already recovered, this gives the uniqueness statement for a solution to 
equation (14. lip . The corresponding stability estimate can be also obtained. 

Finally, let us discuss the case of constant material parameters A, fi, p, vi, . . . , 1/4. In 
this case, to solve equation (14. lip , we do not need to measure ray integrals (L/)(7, 77) 
for all lines 7 of M 3 . Indeed, as is proved in |LSj . three families of lines are sufficient to 
recover the trace free part F from the data KF, each family consists of all lines parallel 
to a coordinate plane. To solve equation f l4.15p . it suffices to know the right-hand side for 
all lines 7 parallel to a plane. 

References 

[CJ] V. Cerveny and J. Jech. Linearized solutions of kinematic problems of seismic body waves in 

inhomogeneous slightly anisotropic media. </. Geophys. 51 (1982), 96-104. 

[HU] S. Hansen and G. Uhlmann. Propagation of polarization in elastodynamics with residual 

stress and travel times. Math. Ann. 326 (2003), 563-587. 

[Ho] A. Hoger. On the determination of residual stress in an elastic body. J. Elasticity 16 (1986), 

303-324. 



TOMOGRAPHY OF SMALL RESIDUAL STRESSES 



17 



[IWY1] V. Isakov, J.N. Wang, and M. Yamamoto. Uniqueness and stability of determining the resid- 
ual stress by one measurement. Comm. Partial Differential Equations 32 (2007), 833-848. 

[IWY2] V. Isakov, J.N. Wang, and M. Yamamoto. An inverse problem for a dynamical Lame's system 
with residual stress. SI AM J. Math. Anal. 39 (2007/08), 1328-1343. 

[IMN] S.A. Ivanov, C.S. Man, and G. Nakamura. Recovery of residual stress in a vertically hetero- 
geneous elastic medium. IMA J. Appl. Math. 70 (2005), 129-146. 

[Kr] Kravtsov Yu.A. "Quasi-isotropic" approximation of geometric optics. Dokl. Acad. Nauk. 

SSSR, 183 (1968), no. 1, 74-76 (in Russian). 

[KO] Kravtsov Y.A. and Orlov Y.I. Geometrical Optics in Inhomogeneous Media. Springer- Verlag, 

Berlin, 1990. 

[LS] W. Lionheart and V. Sharafutdinov. Reconstruction algorithm for the linearized polarization 

tomography problem with incomplete data, in Imaging Micro structures: Mathematical and 
Computational Challenges, Ed. Habib Ammari and Hyeonbae Kang, Contemporary Mathe- 
matics 494 (2009), 137-160. 

[Ma] C.-S. Man. Hartig's law and linear elasticity with initial stress. Inverse Problems 14 (1998), 

313-319. 

[Mu] R. Mukhomctov. On a problem of reconstructing Riemannian metrics. Siberian Math. J 22 

(1982), 420-433. 

[Ra] L. Rachele. Uniqueness in inverse problems for elastic media with residual stress. Comm. 

Partial Differential Equations 28 (2003), 1787-1806. 

[Rol] R.L. Robertson. Boundary identifiability of residual stress via the Dirichlet to Neumann map. 

Inverse Problems 13 (1997), 1107-1119. 

[Ro2] R.L. Robertson. Determining residual stress from boundary measurements: a linearized ap- 

proach. J. Elasticity 52 (1998), 63-73. 

[Sh] V. Sharafutdinov. Integral Geometry of Tensor Fields. VSP, Utrecht, the Netherlands, 1994. 

[Shi] V. Sharafutdinov. Ray transform and some rigidity problems for Riemannian metrics. The 

IMA Volumes in Mathematics and its Applications, 137 (2003) Geometric Methods in Inverse 
Problems and PDE Control, eds. C. Croke, I. Lasieska, G. Uhlmann, M. Vogelius, Springer, 
215-238. 

[Sh2] V. Sharafutdinov. Slice-by-slice reconstruction algorithm for vector tomography with incom- 

plete data. Inverse Problems 23 (2007), 2603-2627. 

Sobolev Institute of Mathematics, Russia 
E-mail address: sriaraf@matri.nsc.ru 



Department of Mathematics, National Taiwan University, Taipei 106, Taiwan 
E-mail address: jnwang@math.ntu.edu.tw 



